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Abstract. Diffusion constants Dr and of ring and linear polymers of 
the same molecular weight in a good solvent, respectively, have been evaluated 
through the Brownian dynamics with hydrodynamic interaction in which no bond 
crossing is possible. The ratio C = Dr/D^, which should be universal in the 
context of the renormalization group, has been estimated as C = 1.14±0.01 for the 
large- ./V limit. It should be consistent with that of synthetic polymers, while it is 
smaller than that of DNAs such as C ~ 1.3. We also perform the same simulation 
through Brownian dynamics with hydrodynamic interaction where bond crossings 
are possible, and obtain almost the same estimate for the ratio C. 
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1. Introduction 

Recently, there has been much progress in experimental techniques associated with ring 
polymers. Ring polymers of large molecular weights are synthesized not only quite 
effectively [T] but also with small dispersions and high purity [H[3]. Diffusion constants 
of linear, relaxed circular and supercoiled DNAs have been measured quite accurately 
[4j. Furthermore, hydrodynamic radius of circular DNA has also been measured [5]. 
The developments are quite remarkable. In fact, it used to be considered quite difficult 
to synthesize ring polymers of large molecular weights. It has now become quite 
interesting to evaluate numerically dynamical or conformational quantities of linear 
and ring polymers that can be measured in experiments. 

It should be nontrivial how linear and ring polymers with the same molecular 
weight in solution may have different dynamical or conformational properties. In fact, 
the excluded volume effect should play a more significant role for ring polymers than 
for linear polymers, since the average distance among monomers is smaller due to 
the constraint of closing two ends [BJ. Moreover, in a dilute solution, the topology 
of a given ring polymer is conserved under thermal fluctuations [7] and represented 
by a knot. Topological constraints may lead to nontrivial statistical mechanical or 
dynamical properties of ring polymers [SI HSl HU H21 HH US ESI H3 - 

In the paper we discuss diffusion constants Dr and Dl of ring and linear polymers 
in good solution, respectively, via the Brownian dynamics with hydrodynamic 
interaction in which bond crossings are effectively prohibited. Here the ring and 
linear polymers have the same molecular weight, and we calculate diffusion constants 
for several different values of the number of segments, N, for 5 < N < 50. We then 
calculate the ratio C — Dr/Dl and compare it with the values measured in some 
experiments and other theoretical values. This gives a test for the validity of dynamical 
models of ring and linear polymers. In fact, it is suggested from the renormalization 
group argument that the large- -/V limit of C should be universal among some class of 
polymer models. Hereafter we call the Brownian dynamics with no bond crossings 
dynamics A. 

Furthermore, we also perform the Brownian dynamics with hydrodynamic 
interaction under almost the same molecular potentials as dynamics A except the 
parameters of the FENE (finite extensible non-linear elongational) potential which 
determine the maximal distance between neighboring monomers. We set the maximal 
distance larger so that bond crossings are allowed. We call it dynamics B. Dynamics B 
has precisely the same potential parameters as that of Ref. [18] . We have found that 
bond crossings occurred for dynamics B, checking the topology of the ring polymer 
by calculating some knot invariants at every time step of the Brownian dynamics. 

Simulation results of both dynamics A and B should be important. In fact, there 
have been several simulation results obtained and accumulated for dynamics B p~8[ 119]. 
We may compare the present simulation with previous ones. In this sense, dynamics 
B is a standard algorithm in the Brownian dynamics. Furthermore, dynamics A is 
important since it preserves the initial topology of a ring polymer. 

The present study should be useful for making explicit connections between 
experimental and theoretical results of dilute solutions of ring polymers. In fact, for 
dilute ring-polymer solutions, even some fundamental properties such as the effects of 
topological constraints have not been clearly confirmed in experiments, yet. Through 
simulations, we can study the effects of topological constraints, which can be checked 
in experiments. 
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The content of the paper consists of the following: In section 2, we briefly explain 
the simulation method. In section 3, we discuss three simulation results. In subsection 
3.1 the ratio of the mean square radii of gyration of ring and linear polymers, g, are 
evaluated numerically. The value of g for dynamics B is consistent with the lattice 
simulation result, while that of dynamics A is larger than the standard one. We 
confirm it also by the Monte-Carlo simulation. It should thus be an interesting future 
problem to evaluate the ratio g for larger values of N. In subsection 3.2, we discuss 
the ratio C both for dynamics A and B. We find that the estimates of C are given by 
almost the same value both for dynamics A and B. Interestingly, the estimate of C is 
consistent with a theoretical value given by a perturbation theory, while it is different 
from that of the renormalization group calculation in one-loop order. However, we 
should note that a one-loop order evaluation could give only a rough estimate and 
multi-loop corrections could improve it. 

2. Simulation method 

The ring polymer molecule is modeled as a cyclic bead-and-spring chain with N beads 
connected by N FENE (finite extensible non-linear elongational) springs with the 
following force law: 



where r — \r\. Let us denote by 6 the unit of distance. Here we assume that the 
average distance between neighboring monomers is approximately given by 6. We set 
constants H and r max as follows: H = 30/csT/6 2 and r max = 1.36 for dynamics A, 
and H = iksT/b 2 and r max — 10b for dynamics B. We assume the Lennard- Jones 
(LJ) potential acting among monomers as follows. 



Here ry is the distance of beads i and j, and 6lj and ctlj denote the minimum energy 
and the zero energy distance, respectively [19] . We set the Lennard- Jones parameters 
as ulj — 0.86 and €lj — O.lfcsT so that they give good solvent conditions [T5]. Here 
ks denotes the Boltzmann constant. 

We employ the predictor-corrector version [20] of the Ermak-McCammon 
algorithm [21 for generating the time evolution of a ring polymer in solution. The 
details are given in Appendix A. The hydrodynamic interaction is taken into account 
through the Ronte-Prager-Yamakawa tensor [551 H3] where the bead friction is given 
by £ = 6Tn] s a with the bead radius a = 0.2576 and a dimensionless hydrodynamic 
interaction parameter h* = ((/6TTr] s )y/H/nkBT = 0.25. 

In the present simulation, physical quantities are given in dimensionless units 
such as in Ref. [19]. We divide length by 6, energy by ksT and time by (Jo 2 /ksT. Let 
us indicate dimensionless quantities by an asterisk as superscript. We have H* = 30, 
r max = 1-3 f° r dynamics A, and H* — 3, r* max — 10 for dynamics B. We take the 
simulation time step At* = 10~ 4 . 

When we evaluate the mean square radius of gyration and the diffusion constant 
for ring and linear polymers through the Brownian dynamics with hydrodynamic 
interaction, we keep each run long enough so that the diffusion constant approaches 
its equilibrium value. For instance, in the case of linear polymers of N = 45 of 
dynamics A, we have performed 9.4 x 10 5 time steps for each run. After the average 



F(r) = -Hr/(l - r 2 /r 2 max ) , 



(1) 




(2) 
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value of the diffusion constant approaches some equilibrium value, we start sampling 
the data and pick up one conformation out of every 18, 800 time steps. Then, the 
diffusion constant evaluated at 740,000th time step is given by 9.256 x 10~ 2 , while 
that at 940,000 th time step is given by 9.276 x 10" 2 . The difference 0.020 x 10~ 2 is 
smaller than their probable error 0.068 x 10~ 2 . 

3. Simulation results 

3.1. Ratio of the mean square radii of gyration 

The mean square radius of gyration (Rq) of a polymer consisting of N monomers is 
defined by 

1 N 

<^> = 4£<(r;-r G ) 2 > 

i=i 

where fj denote the position vectors of monomers for j = 1,2, ...,N, and re the 
position vector of the center of mass of the polymer. The symbol (A) denotes the 
ensemble average of physical quantity A. 

Let us discuss the estimates of the mean square radius of gyration for ring 
and linear polymers, (R g )r and (Rq)l, respectively, obtained by dynamics A and 
dynamics B. They are plotted in Figure 1 against the number of segments N in the 
double logarithmic scales. It is clear that they are fitted well by straight lines. It 
seems that the iV-dependence is close to that of the asymptotic behavior, although 
the number of segments N are not very large, yet. As we shall discuss later, it 
is probably due to the effect of the off-lattice molecular potentials employed in the 
dynamics. Thus, as a fitting formula, we employ the large-iV asymptotic behavior 
of the mean square radius of gyration: (Rq) — AN 2v . The estimates of the fitting 
parameters, An and vr for ring polymers, and Al and i/l for linear polymers, are 
given in the caption of Figure 1. 

Let us now define the geometric shrinking factor g by |24j 

g = (R g )r/(Rq) l . (3) 

We assume that exponent v should be the same for ring and linear chains, i.e. vr = vl- 
We thus have the following fitting formula with three parameters: 

g = 9oc (l + B g N- A *). (4) 

Applying we have g^ = 0.559 ± 0.007 for dynamics A and g^ = 0.535 ± 0.002 for 
dynamics B, as shown in Figure 2. 

The estimate of g value for dynamics B, — 0.535 ±0.002, should be consistent 
with the Monte Carlo simulation using the bond fluctuation model [25] . Interestingly, 
however, the estimate of goo for dynamics A is larger than that of dynamics B even if 
we take into account their errors. The enhancement of value goo in dynamics A should 
be due to the potential forces. In fact, we have confirmed that almost the same value 
of goo is obtained by the Monte-Carlo simulation with the same molecular potentials 
as dynamics A. Therefore, we conclude that it is due to the potential forces employed 
in dynamics A. Here, the potential function of the Monte-Carlo simulation of linear 
chains is given by the following. 

AT-l N 

- MHr 2 max ln[l - (r M+1 /w) 2 ] + 4e LJ ^T^j/r^) 12 - (a LJ /r y -) 6 ] . 
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N 






g = {H z g )r/(H z g)l 


6 


0.413 ±0.001 


0.705 ±0.003 


0.586 ±0.003 


8 


0.606 ±0.001 


1.046 ±0.004 


0.579 ±0.003 


11 


0.916 ±0.002 


1.605 ±0.007 


0.571 ±0.004 


15 


1.364 ±0.002 


2.404 ±0.012 


0.567 ±0.004 


20 


1.961 ±0.005 


3.431 ±0.018 


0.571 ±0.004 


27 


2.842 ±0.008 


5.082 ± 0.024 


0.559 ±0.004 


36 


4.028 ±0.011 


7.182 ±0.037 


0.561 ±0.004 


45 


5.260 ±0.016 


9.335 ±0.053 


0.563 ±0.005 



Table 1. Dynamics A (no bond crossing). Mean square radii of gyration for 
linear and ring polymers, (Rq) r and (R^f)L, and the g values. Applying the 
least square method for {R? g )l = A l N 2u l and (Rq)r = A R N 2v R, respectively, 
the following estimates are obtained: 2u L = 1.288 ± 0.002, Al = 0.072 ± 0.001; 
2v R = 1.270 ± 0.001, A R = 0.043 ± 0.001. In all the tables, errors are given by 
probable errors. 



For ring chains, we add a term of rjv,i due to the periodicity. 

Here we note that we have employed the symbol g^ for the fitting parameter, 
expecting that it should suggest the asymptotic value of g. However, in order to 
evaluate the true asymptotic value of g 7 we have to perform simulations for larger 
values of N. It should thus be an interesting future problem whether the enhancement 
of value g should be relevant to the asymptotic value of g or not. 

According to the one-loop renormalization group calculation [26] g is given by 

5oo = cxp(13/96)/2 = 0.573. (5) 

The value ([5|) is larger than the estimates, g^ = 0.559 ± 0.007 for dynamics A and 
5oo = 0.535±0.002 for dynamics B. Thus, the one-loop calculation does not explain the 
estimate of g^ for dynamics A or B. However, we should note that it is possible that 
the one-loop RG result gives only a crude approximation, and higher-order calculation 
improves the g value. Here we note that for the e-expansion of the n-vector model, 
higher-order terms have been calculated in order to evaluate universal quantities |27j . 
Thus, the multi-loop corrections should be important, although the one- loop correction 
|26j is based on Fixman's cluster expansion [28 and it is not clear whether one can 
extend it. 

Through perturbation calculation, g was estimated in terms of the excluded- 
volume parameter z as follows [29l [30] : 



1 

g= 2 



. 7T 134 \ 

1 + z + 

'2 105 



(6) 



The value of g is dependent on the excluded-volume parameter, z. In order to have 
g ks 0.53, we have z ~ 0.20. Here we note that the value z depends on the number 
of segments N. In order to have z « 0.20 we have to adjust many model parameters. 
Thus, it should be practically impossible to give good estimates of g by making use 
of the perturbation theory. 

3. 2. Ratio of diffusion constants 

Let us recall that the diffusion constant of a polymer is defined by the following: 

D= lim l((rb(t)-rb(0)) 2 ) (7) 
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Figure 1. Mean square radius of gyration (-Rq)l and {Rq)r for dynamics A 
( no bond crossing, nbc) and B (allowed bond crossings, be). For dynamics A, 
data, noints of linear nolvmers are shown bv A. where rinff nolvmers bv H. where 
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Figure 2. The ratio g = (R 2 G ) R / (R 2 G ) L versus N with the fitting curve Q. For 
dynamics A shown by ♦(nbc), we have = 0.559 ± 0.007, B g = 0.402 ± 0.438 
and A„ = 1.173 + 0.706. Here v 2 = .3.3 for S data noints. For dvnamics R shown 
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Figure 3. Diffusion constants of ring and linear polymers for dynamics A (nbc) 
depicted by Aand I, respectively, for dynamics B (be) ♦and V, respectively. The 
horizontal axis denotes the number of segments. N. 
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g = (H g )b/(B; g )l 


6 


0.868 ±0.002 


1.505 ±0.005 


0.577 ±0.003 


8 


1.221 ±0.003 


2.151 ± 0.008 


0.568 ±0.003 


11 


1.774 ±0.004 


3.149 ±0.011 


0.563 ±0.003 


15 


2.550 ±0.005 


4.583 ±0.016 


0.556 ±0.003 


20 


3.587 ±0.008 


6.404 ±0.022 


0.560 ±0.003 


27 


5.172 ±0.012 


9.389 ±0.037 


0.551 ±0.003 


36 


7.155 ±0.017 


12.821 ±0.054 


0.558 ±0.004 


45 


9.306 ±0.023 


17.264 ±0.070 


0.539 ±0.004 



Table 2. Dynamics B (allowed bond crossings). Mean square radii of gyration 
for linear and ring polymers, (Rq)r and {Rq)l, and the g values. Applying the 
least square method for (Rq)l = A l N 2u l and (Rq)r = A R N 2v R, respectively, 
the following estimates are obtained: 2u L = 1.204 ± 0.002, At = 0.175 ± 0.001; 
2v R = 1.179 ± 0.001, A R = 0.105 ± 0.001. 



Here fc(t) denote the position vector of the center of mass of the polymer. Making 
use of ([7]) we have evaluated the diffusion constant of ring and linear polymers, Dr 
and Dl, respectively, through dynamics A and B. 

According to the Einstein relation, the diffusion constant of a polymer should 
be given by D = ksT/C, where C is given by Q = Gitr/Ru with viscosity rj and the 
hydrodynamic radius Rh- Let us assume that the hydrodynamic radius Rh has the 
same asymptotic scaling behavior with the square root of the mean square radius of 
gyration: sj {R 2 G ) oc N v . Thus, in a dilute solution, we have the following large-iV 
behavior: 

D = ksT oc N~ v . (8) 
GnnRH 

Taking the analogy of the large- N behavior JHJ), we introduce the following fitting 
formulas: Dr = A{Dr)N- u{ - d ^ and D L = A{D l )N- u{Dl \ Applying them to the 
data of Table 2, we have the estimates as shown in the caption of Figure 3. The fitting 
curves are shown in Figure 3. The estimates of v(Dr) and v(Dif) are consistent with 
the expected iV-dependence: Dr, Dl oc N~ v with v w 0.59. 

Thus, formula (jHJ) gives good fitting curves to the graphs of the diffusion constants 
Dr and D^ versus N, and the estimates of the exponents v(Dr) and v(Dl) are at 
least roughly in agreement with the SAW exponent v^asn = 0.588, although the large- 
./V behavior ||HJ) should be valid only when N is asymptotically large enough. It is 
likely that N — 50 is not large enough to investigate any asymptotic behavior of the 
diffusion constants. 

It is clear from Tables 3 and 4 that the estimates of C are almost the same for 
dynamics A and B. Here, in Figure 4, the C values are plotted against the number of 
segments N for dynamics A and B, respectively. We also observe that the estimates 
of C are independent of the number of segments, N. In fact, it is also the case with 
the experimental results of DNAs [1] . 

Let us assume again that exponent v should be the same for the diffusion constants 
of ring and linear chains, Dr and Dl, respectively. Applying the fitting formula 

C = Coo (1 + B C N- Ac ) , (9) 

we obtain the following estimate: — 1.14±0.01 for dynamics A; Coo = 1.11±0.01 
for dynamics B. The fitting curves are shown in Figure 4. 
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N 


Dr 


D L 


C = D R /D L 


6 


0.404 ±0.005 


0.368 ±0.005 


1.099 ±0.028 


8 


0.350 ±0.004 


0.308 ± 0.004 


1.137 ±0.029 


11 


0.284 ± 0.004 


0.258 ±0.003 


1.098 ±0.028 


15 


0.244 ±0.003 


0.218 ±0.003 


1.123 ±0.028 


20 


0.199 ±0.003 


0.182 ±0.002 


1.095 ±0.027 


27 


0.172 ±0.002 


0.152 ±0.002 


1.120 ±0.028 


36 


0.142 ±0.002 


0.133 ±0.002 


1.071 ±0.024 


45 


0.131 ±0.002 


0.120 ±0.001 


1.087 ±0.026 



Table 3. Dynamics A (no bond crossing): Diffusion constants of ring and linear 
polymers, D R and Dl, and the C values. Each estimate is derived from the 
average over more than 2,000 runs. We have A(D R ) = 1.011 ± 0.012, u(D R ) = 
0.601 ± 0.004; A{D L ) = 0.938 ± 0.107, v{D L ) = 0.610 ± 0.004. 



N 


D R 


D L 


C = Dr/ Dl 


6 


0.341 ±0.003 


0.316 ±0.002 


1.078 ±0.016 


8 


0.292 ± 0.002 


0.265 ±0.002 


1.102 ±0.016 


11 


0.236 ±0.002 


0.214 ±0.002 


1.104 ±0.016 


15 


0.206 ±0.001 


0.180 ±0.001 


1.140 ±0.016 


20 


0.166 ±0.001 


0.150 ±0.001 


1.109 ±0.016 


27 


0.139 ±0.001 


0.127 ±0.001 


1.100 ±0.017 


36 


0.117 ±0.001 


0.104 ±0.001 


1.122 ±0.018 


45 


0.102 ±0.001 


0.093 ±0.001 


1.107 ±0.017 



Table 4. Dynamics B (allowed bond crossings): Diffusion constants of ring and 
linear polymers, D R and D^, and the C values. Each estimate is given by the 
average over more than 4,000 runs. We have A(D R ) = 1.138 ± 0.022, v(D R ) = 
0.575 ± 0.007; A(D L ) = 1.138 ± 0.022, v(D L ) = 0.610 ± 0.004. 



According to the one-loop renormalization group calculation in the presence of 
both hydrodynamic and self-avoiding interactions |31[ 132] , a universal ratio C is given 

by 

Coo = lim Dr/D l = exp(3/8) = 1.454. (10) 

N— >oo 

The value (fT0|) is much larger than the estimate of C = 1.14±0.01 for dynamics A and 
C = 1.11±0.01 for dynamics B. As in the case of the g value, it is possible that the one- 
loop order result gives only a crude result. Thus, higher-order RG corrections should 
be important. Here we note that the one-loop calculation was performed through the 
conformation-space renormalization-group approach [28j . and it would be nontrivial 
to calculate higher order corrections. 

Some years ago, the ratio C has been estimated by the perturbative calculation 
in terms of the excluded- volume parameter z |33j : 

The value of C is rather constant with respect to z. We have 1.178 at z = 0, and 
1.165 at z = oo. It is interesting to note that the theoretical value (flTj) is rather close 
to the simulation value, C — 1.14 ± 0.01. Thus, the perturbative calculation gives 
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a theoretical value consistent with the simulation result although the validity of the 
perturbation theory is not clear. 

Diffusion constants Dr and have been measured in several experiments. We 
observe a tendency that for synthetic polymers C is given by 1.1 to 1.2, while for linear 
and circular DNAs it is roughly given by 1.3. For instance, it is estimated for relaxed 
circular DNAs as C = 1.32 ± 0.014 [4j. For synthetic polymers through scattering 
experiments, C = 1.1 ~ 1.2 [34] and C = 1.07 ~ 1.15 [35]. Here we note that in Ref. 
[36] C is estimated as a little larger value than in other synthetic polymer experiments. 

We thus conclude that the present model of ring and linear polymers should be 
valid for synthetic polymers, while for relaxed circular DNAs some additional potential 
energy might be important. 

4. Conclusion 

In the present model of the Brownian dynamics both for dynamics A and B, the 
estimate of C = Dr/Dl should be consistent with that of synthetic polymer 
experiments, while it is smaller than that of DNA experiments. The difference in 
the ratio C between synthetic polymers and DNAs may be due to some additional 
potential functions arising from the closed DNA double strands. 
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Figure 4. Ratio C = Dr/Dl versus TV. For dynamics A (4) with fitting curve 
(5), Coo = 1.14±0.01, -B c = 0.094±76.4 and A c = 2.98±454.97. Here x 2 = 6.0 
for 8 data points. For dynamics B (v), C*oo = 1.11 ± 0.01, B c = 161.8 ± 2010.3 
and A c = 4.78 ± 6.94. Here X 2 = 5 for 8 data points. 
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Appendix A. Algorithm of the Brownian dynamics 

In the paper we have simulated linear and ring polymers in a good solvent with 
hydrodynamic interaction by the revised version of the Brownian dynamics [21] with 
respect to the first-order predictor-corrector [20] , 

Let us explain the original version of the Brownian dynamics [21] . We consider N 
Brownian particles in a solvent of temperature T with hydrodynamic interaction. The 
position of the ith Brownian particle, fi, at time t + At is calculated by the following 
equation: 

Ar- = n{t + At) - f t (t) = £ ^ + £ + R t (At) (A.l) 

3 i j B 

for i,j — 1, 2, N. Here, D.y denote the diffusion tensor, Fj the force acting on the 
jth particle, which we shall specify shortly. Ri(At) denote random numbers obeying 
the Gaussian distribution with (Ri(At)) = and (Ri a (At)Rjp(At)) = 2Dij5 a f3At. 

We derive (|A.1[) as follows. First, we consider the Fokker-Planck equation of N 
Brownian particles in a solvent 

¥ =EEfeA^-^iW) (A.2) 

i j j 

where W = W(fi, ..,rjsr,t) is the distribution function for the configuration space of 
the N particles. We can show that the distribution function is given by the multi- 
variable Gaussian distribution if the initial configuration of the N Brownian particles 
is given by W(r\, r^, 0) = Jli ^('^ ~ T 1')- Up to the first order of At, the average 
value and the variance-covariance of the Gaussian distribution, respectively, are given 
by the following: 

(An) =E l^-^*), (A.3) 

(Ar la Ar JI3 ) = 2D l:j S a0 At . (A.4) 

We thus obtain equation (|A. 1|) from the conditions that the difference of the position 
vector Ari = fi(t + At) — fi (t) should satisfy the average value (|A.3[) and the variance- 
covariance (IA.4[) . Here we remark that we can obtain the same average value (|A.3[) and 
the variance-covariance (|A.4[) by integrating the Langevin equations of N Brownian 
particles. 
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Let us now formulate the diffusion tensor and the force acting on the Brownian 
particles. We employ the Ronte-Prager-Yamakawa tensor as the diffusion tensor 



D H = 7r—-h ( for * = 3) ( A - 5 ) 



kT 
6ir(a 
kT 

Da = 

3 8ir( n 



{E+ tM$L\ + ^(\E-- 



ij ij ij 



(iovi^j) (A.6) 



Here a denotes the radius of a bead and C the hydrodynamic friction. For the force, 
we assume the Lennard- Jones force and the FENE spring force. The Lennard- Jones 
potential is given by 

^4 £LJ ((^) 12 -(^) 6 ) (A.7) 

where r is the distance between two particles, ctlj the zero-energy distance and clj 
the energy at distance a. We give ctlj = 0.86 and €lj = 0-lkgT for simulation in a 
good solvent. The potential of the FENE spring force is given by 

Vfene = -\rl iax H\n[l - (-^-) 2 ] (A.8) 

^ 'max 

where r is the distance between a pair of neighboring particles, H the spring constant 
and r max the maximal distance between neighboring particles. For dynamics B, we set 
H* = 3.0 and r* nax = 10.0, which are given in [Tj|]. For dynamics A, we set H* = 30.0 
and r* max — 1.3, as shown in Ref. [37]. In this model no bond crossing should be 
possible due to the strong spring constant and the small maximal distance between 
neighboring particles. Here we note that dimensionless parameters and variables are 
obtained by dividing length, time and energy by b, Cfb 2 jkT and kT, respectively. 

The first-order predictor-corrector version [20] of the Ermak and McCammon 
algorithm [3T] is given as follows. When initial positions of all particles are given, 
we calculate the positions at the next time step as follows. First, we calculate the 
diffusion tensor and the force, i.e. £)°. and F®, respectively, making use of (|A.6[) . 

(|A.7[) and (|A.8[) . Second, we calculate the positions of all particles, r i , by (|A.1|) with 
respect to D° and Ff. Third, using r { , we again calculate the diffusion tensor and 

the force, and denote them by and F i , respectively. Finally, we calculate the 
position of the ith particle at the next time step as follows. 

An = n{t + At) - r?(t) = AtJ2^D^ + ^Ld' 13 ) 

j j j 

+ AtJ2 \ {P%F* + Dyfy/ksT + Rj (A.9) 

3 

Here Rj obey the Gaussian distribution where the average value is zero and the 
variance-covariance is given by the following: 

(R ia R jf3 ) - 2[^{D%+D tj )]5 a0 At. (A.10) 



